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Abstract —  This  paper  presents  a  general  method  for 
investigating  the  unsteady  aerodynamics  of  flapping  wings  for 
micro  air  vehicle  application.  For  this  purpose,  a  dynamically 
scaled  robotic  flapper  was  designed  and  fabricated  which 
can  flap  the  wings  in  a  desired  kinematic  pattern.  A  quasi¬ 
steady  aerodynamic  model  and  wing  testing  methodology 
was  developed  based  on  unsteady  aerodynamic  mechanisms. 
This  model  additionally  accounts  for  the  wing  twisting.  The 
experimental  results  show  a  good  agreement  with  published 
data.  24  kinematic  patterns  were  tested  and  the  quasi-steady 
aerodynamic  model  compares  well  with  the  experimental 
results.  The  focus  of  the  present  work  is  on  hovering  flight, 
however,  the  methodology  is  general  and  can  be  extended  to 
slow  forward  flight  in  future. 

I.  Introduction 
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The  creation  of  flapping  wing  micro  air  vehicles  (FW- 
MAV)  is  a  challenging  problem.  Flapping  wing  flight 
offers  high  maneuverablity  and  the  capability  to  hover  as 
witnessed  in  insects  and  hummingbirds.  These  properties 
make  FWMAVs  more  suitable  for  micro  air  vehicle  mis¬ 
sions  such  as  reconnaissance  and  surveillance,  specially, 
in  confined  areas.  The  aerodynamics  of  flapping  wings, 
such  as  that  of  insects  and  hummingbirds,  is  unsteady. 
The  flow  over  the  wings  is  a  function  of  time  and  this 
makes  unsteady  aerodynamics  extremely  complex.  Unlike 
conventional  fixed  and  rotary  wing  vehicles,  the  flapping 
wing  aerodynamics  is  still  a  largely  unexplored  area. 

As  shown  in  Fig.  1,  the  aerodynamic  module  is  funda¬ 
mental  to  the  design  process  of  a  FWMAV.  The  module 
takes  the  wing  and  body  kinematics  as  inputs  and  gives 
the  aerodynamic  forces  and  torques.  These  are  then  used 
to  compute  the  rigid  body  dynamics,  navigation  and  control 
algorithms,  and  to  perform  design  optimization. 

Due  to  the  complexity  of  solving  Navier- Stokes  equation 
[2]  for  flow  around  flapping  wings  and  possible  inaccu¬ 
racies  in  the  theoretical  modeling,  we  have  selected  the 
experimental  method  to  determine  the  aerodynamic  forces 
and  moments  based  on  blade  element  analysis.  Experi¬ 
mental  investigation  of  flapping  wing  aerodynamics  based 
on  fruit  fly  (Drosophila)  kinematics  has  been  reported 
[9].  Flow  visualization  experiments  using  scaled  hawkmoth 
wings  were  performed  [7,  8].  These  experiments  led  to 
the  discovery  of  certain  unsteady  aerodynamic  mechanisms 
that  are  responsible  for  high  lift  in  biological  flying  species. 


Fig.  1.  The  architecture  of  FW-MAV  Design 


In  our  work,  we  conduct  experiments  on  flapping  wings 
using  a  robotic  flapper.  However,  our  focus  is  the  FWMAV 
aerodynamics  and  design.  Therefore,  we  keep  the  wing 
kinematics  to  be  very  general  and  use  a  generic  insect-like 
wing  for  testing.  We  have  taken  into  account  the  effects  of 
wing  twist  along  the  span,  since  wings  of  large  insects  and 
hummingbirds  show  marked  twisting  compared  to  small 
insects  such  as  Drosophila.  The  moments  and  the  location 
of  force  on  the  wing  are  also  determined  experimentally 
using  a  six-axis  force  torque  sensor.  We  also  present  a 
mathematical  model  of  flapping  wing  aerodynamics  which 
constitutes  the  aerodynamic  module. 

The  main  goal  of  this  paper  is  to  present  a  method 
for  determining  the  force  coefficients  to  be  used  in  the 
aerodynamic  model.  Using  this  methodology,  the  coeffi¬ 
cients  of  a  number  of  wing  shapes  and  geometries  can 
be  catalogued  analogous  to  NACA  airfoil  sections.  This 
information  can  then  be  utilized  for  designing  insect-like 
MAV’s  or  for  comparison  of  aerodynamic  characteristics 
of  different  wing  planforms. 

The  organization  of  this  paper  is  as  follows:  Section  II 
describes  the  flapping  wing  kinematics.  Section  III  presents 
a  derivation  of  the  aerodynamic  model  along  with  the 
assumptions.  Section  IV  outlines  the  experimental  setup. 
In  Section  V,  we  present  a  modification  to  the  aerodynamic 
model  based  on  the  initial  experimental  data.  Section  VI 


explains  the  methodology  of  determining  the  coefficients  in 
the  aerodynamic  model.  In  Section  VII,  the  experimental 
results  are  compared  with  published  data  and  in  Section 
VIII,  we  compare  the  quasi- steady  aerodynamic  model  with 
the  experimental  results.  Finally,  conclusions  are  presented 
in  Section  IX. 

II.  Kinematics 
A.  Wing  motion  analysis 

We  consider  the  left  wing  of  a  typical  insect.  For  the 
right  wing  motion,  a  similar  but  different  rotation  sequence 
is  required  and  will  not  be  considered  here.  As  shown  in 
Fig.  2,  the  wing  frame  (xw,yw,  zw)  can  be  described  by 
three  succesive  rotations  with  respect  to  frame  (xs,ys,  zs ) 
attached  at  the  wing  base.  First,  rotation  about  axis  by 
an  angle  </>,  next  a  rotation  about  the  current  x'  axis  by  an 
angle  0  and  finally  rotation  about  the  current  y"  axis  by 
an  angle  ip .  Therefore,  the  wing  position  is  given  by  the 
body  sequence  3-1-2  rotation  and  by  angles  (</>,  ip).  The 
yw  axis  attached  to  the  wing  is  the  feather  axis. 


where  u  is  the  flapping  speed  in  radians/sec,  <f>n,  0n, 
are  the  flapping  amplitudes  and  (en  and  C^n  are  the 
phase  differences. 

B.  Terminology 
Stroke  plane 

The  plane  defined  by  (xs,  ys )  axes  represents  the  stroke 
plane  as  shown  in  Fig.  2.  The  frame  (xs,  ys,zs)  is  refered 
to  as  the  stroke  plane  frame.  The  motion  of  the  feather 
axis  is  not  necessarily  confined  to  the  stroke  plane  due  to 
6 ,  the  elevation  or  out  of  stroke  plane  angle. 

Planar  flapping 

If  the  out  of  stroke  plane  angle  6  is  zero,  the  motion 
is  called  planar  flapping.  In  this  paper,  we  will  consider 
planar  flapping. 


Fig.  3.  Motion  terminology:  (A)  Rotational  Motion  of  a  wing 
section  about  the  feather  axis  given  by  ip.  (B)  Translational 
motion  of  the  section  about  the  flapping  axis  given  by  cp.  The 
difference  between  flapping  translation  and  linear  translation 
is  shown  in  (C)  and  (D). 
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Fig.  2.  Figure  showing  the  body  sequence  3-1-2  which  gives 
the  three  rotation  angles  (0, 0 ,  ip) 

Wing  angles  can  be  a  represented  by  Fourier  series  given 
below  [10]: 

oo 

4>(t)  =  4>0  +  y2®nSm(Tuot+ (1) 

n= 1 

oo 

9(t)  =  0O +  y2®nSm(riL0t  +  (en),  (2) 

n= 1 

oo 

ip(t)  = -ipo +  y2^nSin(ncot  + J,  (3) 

n= 1 


Linear  translation  [11]  is  the  wing  motion  where  the 
entire  wing  moves  forward  perpendicular  to  the  feather  axis 
in  such  a  way  that  the  velocity  of  the  base  and  the  tip  of  the 
wing  are  the  same  as  shown  in  Fig.  3(D).  Therefore,  every 
section  of  the  wing  moves  along  a  straight  line  path  with 
the  same  velocity.  A  conventional  aircraft  wing  undergoes 
linear  translation.  In  flapping  translation  [11],  the  wing 
tip  rotates  around  the  flapping  axis  fixed  at  the  base  as 
shown  in  Fig.  3(C).  Therefore,  the  velocity  along  the  wing 
increases  linearly  from  the  root  to  the  tip  such  that  a  section 
of  the  wing  translates  along  a  circular  path  as  shown  in  Fig. 
3(B).  Insect  wings  undergo  flapping  translation.  A  section 
of  helicopter  rotor  blade  undergoes  flapping  translation 
since  the  section  follows  a  circular  path.  However,  unlike 
insect  flapping  motion,  it  can  trace  a  complete  circle.  For 
a  wing  section  located  at  a  distance  r  from  the  flapping 
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axis,  the  translational  speed  \ut(tR)\  is  given  by  <pr  as 
shown  in  Fig  4(C)  and  the  translational  acceleration  along 
the  circular  path  is  cpr. 

Flap  amplitude 

For  a  single  frequency  motion,  the  flap  amplitude  is  an 
important  parameter  of  the  wing  kinematics,  denoted  by 
as  shown  in  Fig.  3(B).  The  total  stroke  amplitude  is  labeled 
as  <l>,  which  is  twice  the  flap  amplitude  T>.  T>o  is  the  mean 
stroke  amplitude.  The  mid  —  stroke  is  that  position  of  the 
wing  where  (j>  =  4>0  as  shown  in  Fig.  3(B).  When  <f>o  7^  0, 
the  motion  is  called  asymetric  flapping.  When  4>0  =  0, 
the  motion  is  refered  to  as  symmetric  flapping. 

Rotational  motion 

Wing  rotation  refers  to  the  rotation  of  the  wing  about 
the  feather  axis  through  an  angle  ip  as  shown  in  Fig.  3(A) 
and  given  by  Eq.  (3).  Flapping  wing  typically  rotates  at 
the  extremes  of  the  stroke.  When  the  wing  goes  from 
upstroke  to  downstroke,  the  rotation  is  called  pronation. 
Similarly,  rotation  between  downstroke  and  upstroke  is 
called  supination.  The  rotational  speed  and  rotational 
acceleration  are  ip  and  ip,  respectively. 

Rotational  amplitude 

The  rotational  amplitude  is  denoted  by  T.  As  shown 
in  Fig  3(A),  \k0  is  the  mean  rotational  amplitude.  For 
symmetric  rotation ,  4/q  =  0,  otherwise  the  motion  is 
refered  to  as  asymetric  rotation. 

Angle  of  attack 

The  angle  of  attack  is  the  angle  between  the  chord  at  a 
location  r  from  the  wing  base  and  the  relative  flow  velocity 
ur(r,t).  Mathematically,  a  is  given  by 

u 

a  =  ±tan(  —  )  ,  —  7r/2  <  a  <  7r/2,  (4) 

Uf 

where  ut  and  un  are  the  components  of  ur(r ,  t)  along  the 
normal  n  and  axial  t  axes  fixed  to  the  airfoil  section  shown 
in  Fig.  4(B).  The  db  in  Eq.  (4)  takes  the  sign  of  ut. 

Mean  angle  of  attack  (a) 

We  define  the  mean  angle  of  attack  a  as  the  angle  of 
attack  of  a  section  located  at  a  distance  ^77  along  the 
span  and  having  the  chord  length  equal  to  the  mean  chord 
c  [3].  r 2  is  the  radius  of  gyration  of  the  wing  area  non- 
dimensionalized  by  the  wing  length  R  [3].  Therefore,  a  is 
given  by  Eq.  (4) 

a  =  iarctan  ( ,  —  it/2  <  a  <tt/2  (5) 

\ut(r2R,t)  ) 


V  Feallier  aids 

or  dotation  axis 


Fig.  4.  (A)  Wing  geometry  of  a  typical  insect.  (B)  shows  the 

definition  of  angle  of  attack  at  a  wing  section.  (C)  The  angle 
of  attack  remains  constant  during  the  translational  motion  of 
a  section.  If  translational  speed  \iiT(r,t)\  =  cpr  is  also  held 
constant,  then  the  flow  over  the  section  is  steady-state. 

III.  Aerodynamic  model 
A.  Quasi-steady  state  analysis 

If  the  flow  velocity  at  a  location  on  the  wing  does  not 
change  with  respect  to  time,  it  is  referred  to  as  steady  state 
flow.  In  order  to  maintain  steady  state  flow  over  an  airfoil 
section,  it  is  required  that  the  section  be  placed  in  a  flow 
of  constant  velocity  and  at  a  fixed  orientation,  i.e,  the 
angle  of  attack  a  with  respect  to  time  be  fixed  in  an  inertial 
frame.  The  steady- state  aerodynamic  force  is  given  by 

F  =  CF{a)  1/2  p  UjA.  (6) 

where  F  is  the  aerodynamic  force  such  as  lift  or  drag, 
is  the  relative  flow  velocity,  A  is  the  wing  area  and  p  is 
the  air  density.  Cp(a)  is  the  coefficient  of  the  aerodynamic 
force.  It  is  a  function  of  the  angle  of  attack  a  of  the  wing. 

The  steady  state  force  equation  for  flapping  wing  can  be 
derived  using  the  blade  element  method  (BEM).  The  force 
on  a  strip  of  wing  at  a  distance  r  from  the  flapping  axis  is 
given  by 

dF  =  Cf{ol)  1/2  p  \uT(r,t)\2c(r)dr ,  (7) 

where  Cpia)  is  the  force  ceofficient  of  an  element  of  the 
wing  which  is  a  function  of  the  local  angle  of  attack  a  and 
c(r)dr  is  the  area  of  the  section  as  shown  in  Fig.  4(A).  Eq. 
(7)  can  be  integrated  for  the  entire  wing  as  follows 

rR 

F  =  1/2  p  /  CF(a)  \uT{r,t)\2c(r)dr,  (8) 

Jo 

The  force  coefficient  can  also  be  defined  for  the  entire  wing 
based  on  the  mean  angle  of  attack  a.  Therefore,  Eq.  (8) 
can  be  written  as 


(9) 


rR 

F  =  CF{a)l/2p  /  \uT(r,t)\2c(r)dr. 

Jo 

In  the  hope  of  finding  approximate  analytical  solutions  to 
the  flapping  wing  aerodynamics,  simplified  models  based 
on  quasi-steady  state  assumption  have  been  developed  [15]. 
According  to  quasi-steady  state  assumption,  the  motion 
kinematics  during  flapping  cycle  is  replaced  by  a  series 
of  static  positions  having  instantaneous  velocity  and  angle 
of  attack  [3].  The  force  is  determined  using  Eq.  (7)  or  (8) 
which  is  not  a  function  of  wing  rotation  and  acceleration. 
It  is  only  a  function  of  the  translational  velocity  ur(r,t) 
of  the  wing.  Therefore,  we  refer  to  quasi-steady  state 
force  as  the  translational  force.  In  this  method,  any  time 
dependence  of  the  aerodynamic  force  arises  from  the  time 
dependence  of  the  kinematics  but  not  that  of  the  fluid 
flow  itself.  Ellington  [3]  used  quasi-steady  analysis  to 
investigate  insect  flight  and  determined  that  the  analysis 
underestimated  the  lift  required  to  support  an  insect  during 
hovering  and  concluded  that  a  substantial  revision  of  the 
quasi-steady  method  is  necessary.  He  further  proposed  that 
the  quasi- steady  state  model  must  include  wing  rotation  in 
addition  to  flapping  translation. 

B.  Quasi-unsteady  state  analysis 

Since  Ellington’s  investigation  [3],  several  researchers 
have  provided  more  data  to  support  the  insufficiency  of 
the  quasi-steady  model  (Ennos,  1989a;  Zanker  and  Gotz, 
1990;  Dudley,  1991).  These  developments  have  spurred  the 
search  for  specific  unsteady  mechanisms  and  mathematical 
models  to  explain  the  aerodynamic  forces  on  insect  wings. 
Dickinson  [9], [12]  used  experimental  investigation  to  de¬ 
termine  the  aerodynamics  of  hovering  fruitfly.  According  to 
Dickinson  [12],  the  total  instantaneous  aerodynamic  force 
on  the  wing  can  be  represented  as  a  sum  of  four  force 
components  given  below; 

F'inst  =  Fa  +  F trans  T~  F rot  T~  F wo  (10) 

where  Finst  is  the  instantaneous  aerodynamic  force,  Fa 
is  the  force  due  to  virtual  mass  effect  and  Ftrans  is  the 
instantaneous  translational  force,  Frot  is  the  rotational  force 
and  Fwc  is  the  force  due  to  wake  capture.  Dickinson 
did  not  provide  a  mathematical  form  for  wake  capture 
effect.  However,  wake  capture  effect  was  identified  from 
the  experimental  data.  Delaurier  [1]  presented  a  theoretical 
model  of  flapping  wings  based  on  circulation  theory  of 
lift.  In  Delaurier’ s  model,  the  wake  capture  effect  is  also 
modeled  in  addition  to  rotational  and  virtual  mass  forces. 
Walker  [13]  also  presented  a  semi-empirical  model  of 
flapping  wings  with  an  alternate  mathematical  form  for 
the  rotational  force.  Walker  applied  this  model  to  the 
fruitfly  wing  and  found  good  comparison  with  both  the 
experimental  results  of  fruitfly  experiment  by  Dickinson  et 


al  and  the  CFD-modeled  forces  on  the  virtual  fruitfly  wing 
(Sun  and  Tang,  2002  [16]). 

Based  on  a  study  of  quasi-unsteady  aerodynamic  models 
presented  by  Dickinson  [12],  Delaurier  [1]  and  Walker  [13], 
we  modify  the  quasi-steady  model  given  by  Eq.  (8)  as 
follows 

Ft  =  Fsteady{4>)  H-  FunstQadyicj)')  ip,  (f),  pf) 5  (11) 

where  Ft  is  the  total  instantaneous  force  on  the  wing, 
F steady  is  the  steady  state  or  translational  force  given  by 
Eq.  (9)  and  Funstea(iy  is  the  unsteady  force  which  is  a 
function  of  wing  rotation  and  acceleration.  In  this  paper, 
we  will  refer  to  this  method  as  quasi-unsteady  since  the 
total  force  Ft  is  implicitly  dependent  on  time.  The  current 
focus  is  to  determine  the  mathematical  form  of  Funsteady. 

C.  Unsteady  effects 

Leading  edge  vortex  (LEV)  Force 

When  a  thin  wing  translates  at  a  high  angle  of  attack 
close  to  the  stall,  the  flow  breaks  up  at  the  leading  edge  and 
rolls  into  a  leading  edge  vortex  (LEV).  The  presence  of  this 
vortex  increases  the  circulation  and,  thereby,  the  lift  force 
significantly.  In  conventional  airplane  wings,  this  effect 
occurs  momentarily  before  stall.  However,  in  insects  such 
as  hawkmoths,  the  flapping  translation  motion  stabilizes  the 
vortex  and  it  remains  attached  to  the  wing  during  the  entire 
stroke  [2],  [7].  LEV  was  shown  to  remains  attached  even  in 
rotating  wings  at  high  angles  of  attack  [5],  [6].  Therefore, 
based  on  this  study,  we  conclude  that  LEV  force  is  not  a 
function  of  wing  rotation  and  acceleration  and  it  can  be 
modeled  by  steady-state  Eq.  (9)  which  can  be  written  as 

rR 

Fiev  =  C1(a)l/2p  /  \uT{r,t)\2c(r)dr  =  CTx.  (12) 
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Equation  (12)  can  be  thought  of  as  a  product  of  a 
function  F\  =  1/2  p  \uT(r,t)\2c(r)dr  which  captures 
the  physics  of  Fsteady  =  Fiev  and  a  coefficient  C%  — 
Ci  (a)  that  adjusts  the  magnitude.  Since  Fiev  is  generated 
during  the  translational  phase  of  wing  motion,  it  can  be 
refered  to  as  translational  force. 

Rotational  Force 

If  the  wing  rotates  about  the  feather  axis  with  an  angular 
rate  tp,  a  rotational  circulation  force  is  generated  [12].  In 
this  paper,  we  use  the  mathematical  form  of  rotational 
force  given  by  Walker  [13].  It  says  that  rotational  force 
can  be  modeled  by  selecting  the  flow  velocity  \uT(r,t)\ 
in  Eq.  (12)  at  a  location  T  along  the  chord  as  shown 
in  the  Fig.  5(A).  The  total  flow  velocity  u(r,t)  can  be 
written  as  a  vector  sum  of  translational  velocity  UT(r,t) 
and  rotational  velocity  ur(v ,  t).  The  magnitude  of  fi#(r,  t) 
is  dependent  on  dQ  and  di  which  are  percentage  distances 
along  the  chord  c(r).  The  parameter  di  is  an  unknown 


constant  while  dQ  is  known  from  wing  geometry  and  can 
vary  along  the  span.  Therefore  dQ  =  d0(r).  The  coefficient 
of  rotational  force  C2  appears  as  the  non-dimensional 
parameter  di  —  dQ(r )  in  the  expression  for  \u(r,t)\  which 
can  be  adjusted  to  scale  the  rotational  force.  Therefore, 

C2(r)  =  di  -  da(r).  (13) 

Note  that  if  the  rotational  axis  lies  ahead  of  the  leading 
edge  then  d0(r)  should  be  taken  as  negative.  Therefore  C2 
can  vary  along  the  span  depending  upon  the  wing  geometry 
unless  dQ(r)  =  0.  The  combined  LEV  and  rotational  force 
is  given  as 

fR 

Flev+rot  =  Ci  1/2  p  /  \u{r,t)\2c{r)dr,  (14) 
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where  the  coefficient  C2  appears  in  the  expression  for 
|  u(r,  t)  |.  If  the  wing  has  only  translation,  i.e,  ^  is  zero,  then 
Eq.  (14)  simplifies  to  Eq.  (12)  and  we  get  the  translational 
force  only.  In  Fig.  5(A),  note  that  the  angle  of  attack  a  due 
to  combined  rotational  and  translational  velocity  is  greater 
than  olt  due  to  translational  velocity  alone. 


Fig.  5.  (A)  The  combined  rotational  and  translational/LEV 

effect  on  a  wing  section  located  at  a  distance  ‘r’  from  the 
flapping  axis.  ur(r,t)  is  the  translational  velocity,  UR(r,t) 
is  the  rotational  velocity  and  iqr,  £)  is  the  vector  sum  of 
translational  and  rotational  velocities.  (B)  Virtual  mass  force 
acting  at  a  section  located  at  a  distance  ‘r’  from  the  flapping 
axis.  Here,  un  (r,  t)  is  the  component  of  acceleration  normal 
to  the  wing  surface  in  the  wing  frame  and  ‘  dnrd  is  the  mass 
of  air  assumed  to  be  contained  in  a  cylinder  of  diameter  c(r) 
and  height  dr. 

Virtual  Mass  force 

As  the  wing  accelerates,  it  moves  along  with  it  a  mass 
of  air,  assumed  to  be  contained  in  a  cylinder  with  diameter 
equal  to  the  chord  [1],  [3].  The  acceleration  of  this  mass 


of  air  shows  up  as  a  virtual  mass  force  (see  Fig.  5(B))  and 
can  be  written  as 

dFvirtual  mass  =  dlTi  Un(r:tf  (15) 

where  iin(r,t)  is  the  rate  of  change  of  normal  velocity 
component  at  the  mid-chord  location  in  the  wing  frame  and 
dm  is  the  mass  of  air  enclosed  in  a  thin  cylinder  of  width 
dr  and  a  diameter  equal  to  the  chord  c  at  a  distance  r  from 
the  flappping  axis.  The  mass  of  air  is  pnc2  /  4.  Therefore, 
Eq.  (15)  can  be  written  as 

dFvirtuai  mass  =  —^-un(Vi  t)c{r)  dr.  (16) 

On  integrating  Eq.  (16)  from  root  to  the  tip  of  the  wing, 
we  get  the  total  force  given  by 

Fyirtual  mass  =  6/3——  /  Un{r  ^  t)c{r)  dr  =  C3F3  (17) 
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The  coefficient  of  virtual  mass  C3  is  included  in  order  to 
adjust  the  magnitude  while  the  function  F%  captures  the 
physics  of  virtual  mass  effect. 

D.  Total  Force 

The  total  aerodynamic  force  including  translational 
(LEV  effect),  rotational  and  virtual  mass  effects  can  be 
written  as 

Ft  —Ci^  J  \u{r,t)\2c{r)dr  + 

£3^7-  /  un(r,t)c(r)2dr.  (18) 
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Here,  the  LEV  and  rotational  forces  are  combined  as 
the  first  term  in  the  above  equation.  C\  and  C3  are  the 
coefficients  of  LEV  and  virtual  mass  force  respectively. 
The  coefficient  of  rotational  force,  C2,  appears  implicitly 

in  u(r,  t). 

E.  Assumptions  in  the  model 

•  Superposition  of  steady  and  unsteady  aerodynamic 
effects  holds. 

•  LEV  force  can  be  modeled  by  the  steady- state  aero¬ 
dynamic  equation. 

•  Chordwise-force  due  to  skin  friction  is  ignored.  This 
is  based  on  the  results  of  revolving  wing  experiments 
[5], [6]. 

•  The  total  force  Ft  acts  normal  to  the  chord  through¬ 
out  the  flapping  cycle,  i.e,  we  assume  that  C\Fi, 
F2(CiC2)  and  C3F3  act  normal  to  the  chord  at  every 
section  of  the  wing. 

•  The  total  force  Ft  acts  at  the  mid-chord  location  at 
every  section  of  the  wing. 


IV.  Experimental  Investigation 


A.  Flow  similarity 

In  order  to  determine  the  coefficients  C\,  C2  and  C3 
in  Eq.  (18),  we  conducted  experimental  investigation.  The 
basis  of  experimental  investigation  is  flow  similarity  which 
ensures  that  the  coefficients  are  similar  for  the  actual  and 
the  experimental  wing.  In  order  to  achieve  flow  similarity, 
the  reduced  frequency  parameter  K  along  with  Reynolds 
number  Re  and  wing  geometry  should  match  for  the  actual 
and  experimental  wing.  The  wing  size,  flapping  speed  and 
fluid  medium  can  be  different.  Therefore, 


Ci  =  C2  =  f  (Re,  K,  geometry)  (19) 

It  is  shown  in  [17]  that  the  virtual  mass  force  depends  on 
Reynolds  number  and  wing  geometry 


C3  =  f(Re ,  geometry)  (20) 

The  flow  Reynolds  number  Re  and  reduced  frequency  K 
for  the  case  of  hovering  flapping  flight  are  given  by  [4,  10] 


Re  = 


8<f>A2/ 


K  = 


(21) 


i/A  ’  24>A 

where  v  is  the  kinematic  viscosity  of  the  fluid  medium,  R 
is  the  wing  length,  <f>  is  the  flapping  amplitude,  /  is  the 
flapping  frequency  in  cycles/sec  and  A  is  the  wing  aspect 
ratio.  Reduced  frequency  AT  is  a  measure  of  unsteadiness 
of  the  flow. 


B.  Robotic  Flapper 

In  order  to  conduct  experimental  investigation,  a  robotic 
flapper  was  designed  and  fabricated  at  the  University  of 
Delaware.  It  is  shown  in  Fig.  6.  The  flapper  is  driven 
by  three  independent  servo  motors  designed  to  give  three 
degrees-of-freedom  flapping  motion,  i.e,  (p,  6 ,  vp. 

Flapper  Kinematics 

In  our  experiments,  we  keep  the  out  of  plane  motion 
6  =  0.  This  simplifies  the  kinematics  but  still  retains  the 
features  of  wing  motion.  The  flapper  coordinate  system 
is  shown  in  Fig.  7(A).  For  the  case  of  hovering  flight, 
the  body  of  FWMAV  is  assumed  to  be  stationary  with 
respect  to  the  earth.  The  body  frame  f0(xo,  2/0  5  io)  is  also 
the  inertial  frame.  The  rotation  matrices  between  the  body 
frame  and  frame  f\(x\,yi,  z\)  and  between  fi  and  wing 
frame  f3  (x3,  y3,  z3)  are  given  by: 


Fig.  6.  Figure  shows  the  robotic  flapper  designed  and 
fabricated  at  University  of  Delaware.  It  is  driven  by  three 
independent  servo  motors  and  can  give  3-DOF  flapping  wing 
motion.  A  six-axis  force  torque  sensor  (Nano  17)  from  ATI 
industrial  automation  is  mounted  at  the  base  of  the  wing.  The 
aerodynamic  forces  and  torques  along  with  wing  position  can 
be  seen  real  time  with  the  help  of  this  apparatus. 


the  directions  (£4, 1/4,  £4)  respectively.  The  rotation  matrix 


between  the  wing  frame  f3  and 

sensor  frame  /4  is  given 

by 

/-1  0 

°\ 

=  0  -1 

0  ,  (23) 
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The  lift,  drag  force  is  given  by 


(24) 


0 

Lift 

Drag 


(25) 


where  Lift  and  Drag  are  along  y\  and  z\  direction 
respectively.  Similarly,  the  vertical  and  horizontal  force  is 
given  by. 


(C^  0  S^\  (  0  0  1\ 

/?,'  =  [  S*  0  -C*  ,  R\  =  -S+  0 

V  0  1  0/  V  s4,  c\,,  ° ) 

(22) 

The  frame  f 4(004, 7/4,  £4)  is  the  force/torque  sensor 
frame.  The  Fx,Fy,Fz  are  the  sensor  read  forces  along 


•  (26) 

where  FyL  is  the  horizontal  force  and  Fv  is  the  vertical 
force  along  the  directions  yQ  and  zQ  respectively.  Due  to 
symmetry,  we  expect  the  force  along  the  ocq  direction  to 


Fig.  7.  (A)  Figure  shows  the  sensor  co-ordinate  frame  /4 

and  the  postive  direction  of  force  and  moment  components 
{Fx,Fy,Fz,Mx,My,Mz).  The  direction  of  lift,  drag,  vertical 
and  horizontal  force  are  also  identified.  (B)  Wing  Planform 
used  in  the  experiment. 

cancel  by  the  two  opposite  wings.  The  angular  velocity  of 
wing  with  respect  to  the  earth  frame  is  given  by 

^3/o  =  <j>z0  +  ipz3.  (27) 

Therefore,  velocity  of  air  at  the  point  ’£  on  the  wing  in 
frame  /4  is  given  by 

u(r,t)  =  —  (prsmipx^  —  (pr  cos  ip  +  C2c(r)ip)y4  + 

C2c(r)<p  sin  ip £4.  (28) 

For  the  estimation  of  total  aerodynamic  force  given  by  Eq. 
(18),  the  velocity  and  normal  acceleration  components  are 
required.  For  the  robotic  flapper,  these  are  given  by 

\u(r,  t)\2  =  (<pr  simp)2  +  (pr  cosip  +  C2c(r)ip)2 ,  (29) 

un(r,  t )  =  —  (pr  cos  ip  —  ip  p  simp  +  0.5c(r)  ip),  (30) 

where  the  velocity  of  air  in  the  £4  direction  is  ignored 
because  it  is  in  the  spanwise  direction  and  does  not 
contribute  to  the  force.  For  estimation  of  a ,  given  by  Eq. 
(5),  the  normal  and  tangential  velocity  components  in  the 
2/4  and  x4  directions  respectively  are  given  by 

un{r2R,t)  =  -  (pr2R  cos  ip  +  C2cip),  (31) 

ut(r2R,t )  =  —  pr2Rsinp.  (32) 


If  we  substitute  \u(r,  t)  |2  and  un{r,  t )  given  by  Eqs.  (29) 
and  (30)  into  Eq.  (18),  we  get  the  total  force  as 

Ft  =  C\F\  +  F2[C\,  C2)  +  C3F3,  (33) 

where 

Fi  =  ^  cp2  f  r2c(r)dr  (34) 
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F2  =  |  [2CiC2ip(p  cos  ip  rc(r)2dr  + 

pR 

Clip2  /  c(rf dr]  (35) 

Jo 

F3  =  ^[—(p  cosip  fQR  rc(r)2dr  +  ipcp  sin  ip  fR  c(r)2  dr 

..  rR 

-0.5 ip  /  c(r)3dr]  (36) 
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where  C\F\  is  the  steady-state  or  LEV  force,  C3F3  is  the 
virtual  mass  force.  The  function  F2  captures  the  rotational 
force  and  it  is  an  implicit  function  of  the  coefficients  C\ 
and  C2.  The  integrals  are  a  function  of  wing  geometry. 

C.  Wing  Design 

Size  and  shape:  The  size  of  the  wing  can  be  determined 
from  Eq.  (21)  based  on  Reynolds  number.  Reynolds  num¬ 
ber  of  hummingbird  ( Lampornis  clemenciae )  is  14,300 
with  a  wing  length  of  85  mm,  aspect  ratio  of  8.2,  total 
stroke  amplitude  of  151  deg  and  a  flapping  frequency  of 
23.3  cycles/sec.  A  flapping  wing  MAV  capable  of  hovering 
flight  will  be  similar  to  a  hummingbird  in  size  and  therefore 
the  anticipated  Re  range  for  FWMAV  design  is  between 
12,000-25,000.  Based  on  this  Re  range,  we  selected  a 
generic  insect-like  wing  shape  for  experiment  shown  in 
Fig.  7(B).  It  has  the  following  scaled  dimensions: 

•  Wing  length  R  =  0.58m, 

•  Aspect  ratio  A  =  5.7677, 

•  r2  =  0.5628  (dimensionless). 

•  Scaled  flapping  frequency  =  0.5  cycles/sec 

The  wing  geometry  is  given  in  Table  I.  Each  element  has 
constant  width,  i.e  dr  =  0.05  meters. 

Wing  fabrication 

The  wing  structure  is  made  of  carbon  rods  representing 
the  veins  in  a  typical  insect  wing  or  feathers  in  humming¬ 
bird’s  wing.  Carbon  rods  radiate  from  the  triangular  wing 
base  made  of  balsa  block.  A  cellophane  membrane  was 
attached  to  the  structure  using  cellophane  tape.  The  entire 
wing  assembly  weighs  just  18  grams. 


TABLE  I 

Wing  Geometry  Used  in  the  Experiment 


section 

spanwise  location  T’(in), 

chord  c(r)(m), 

d0{r)% 

1 

0.105 

0.270 

0 

2 

0.155 

0.287 

0 

3 

0.205 

0.284 

0 

4 

0.255 

0.271 

0 

5 

0.305 

0.261 

0 

6 

0.355 

0.254 

0 

7 

0.405 

0.236 

0 

8 

0.455 

0.212 

0 

9 

0.505 

0.180 

-0.12 

10 

0.555 

0.115 

-0.27 

D.  Measurement  of  Force  and  Moment 

We  used  ATI  Industrial  Automation  multi-axis  force- 
torque  sensor  which  can  measure  three  forces  Fx.  Fy.  Fz 
and  three  moment  components  Mx,  My,  Mz.  The  sensor 
can  measure  a  maximum  of  ±12  N  in  the  Fx  and  Fy 
directions  and  ±17  N  in  the  Fz  direction.  The  maximum 
moment  range  in  all  directions  is  ±120  N-mm.  The  reso¬ 
lution  for  all  three  force  components  is  1/1280  N  and  for 
all  three  moment  components  it  is  1/256  N-mm. 

In  order  to  reduce  noise  from  the  data,  we  used  a  simple 
low  pass  digital  filter  [14].  It  is  given  by  the  following 
difference  equation 


y(n)  =  a  y(n  —  1)  ±  (1  —  a)  x(ri),  (37) 

where  x(n)  represents  the  discrete-time  observed  signal 
with  n  =  0,1,  2. .N  at  the  sample  points.  y(n )  is  the 
smoothed  output  and  y{— 1)  =  0.  The  parameter  a  is  a 
weighting  factor  (0  <  a  <  1)  selected  between  0.7  and 
0.8,  depending  on  the  noise  in  the  data. 

V.  Initial  observations 

The  original  force  and  torque  data  along  with  filtered 
data,  in  the  sensor  coordinate  frame  for  one  particular 
kinematic  pattern  is  shown  in  Fig.  8.  This  data  includes 
the  forces  and  moments  due  to  gravity  as  well  as  inertia. 
In  order  to  filter  out  the  aerodynamic  force,  we  used  an 
identical  wing  but  without  the  membrane  (wing  B).  The 
wing  with  the  membrane  is  refered  to  as  wing  A.  The 
inertia  and  gravity  loads  from  wing  B  were  subtracted  from 
wing  A  to  get  the  aerodynamic  force  data. 

Fig.  9  shows  the  aerodynamic  and  inertia  loads  from 
wing  A,  the  inertia  loads  from  wing  B  and  the  aerodynamic 
force  in  the  sensor  Fx,  Fy  and  Fz  directions  for  one  stroke. 
A  ten  degree  polynomial  fit  was  done  on  the  aerodynamic 
data  to  further  smooth  out  the  noise. 

We  know  from  revolving  wing  experiments  [5],  [6]  that 
at  high  angles  of  attack  (10°  and  up),  the  aerodynamic  force 
is  roughly  normal  to  the  wing  surface  and  this  was  the 
assumption  made  in  our  aerodynamic  model.  This  implies 
that  the  sensor  should  only  detect  the  Fy  component  of 
force.  However,  initial  results  revealed  the  presence  of  Fx 
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Fig.  8.  Force  Data  in  sensor  frame  of  wing  A  and  wing  B 
(before  and  after  filtering).  The  unit  of  force  is  ‘N’  (Newtons) 
and  the  unit  of  moment  is  ‘N-mm’ 


and  Fz  components  of  the  force  besides  the  Fy  component 
as  shown  in  Fig.  8.  From  observations  of  the  wing  during 
the  experiment,  we  concluded  that  the  wing  deforms  due  to 
the  aerodynamic  loads.  Both  span  wise  bending  and  twist 
were  present  which  are  discussed  separately  below. 

Spanwise  Bending 

Initially,  the  Fz  data  appeared  simply  as  noise  as  shown 
in  Fig.  8.  However,  after  processing  the  data  for  one 
complete  cycle,  we  observed  that  there  was  a  definite 
aerodynamic  force  in  the  Fz  direction.  Fig.  9  shows  the 
aerodynamic  force  Fz  along  the  sensor  direction  which 
is  along  the  leading  edge  or  feather  axis  of  the  wing  (see 
Fig.  7(A)).  If  we  examine  the  inertia  loads  from  wing  B, 
we  see  that  there  is  a  positive  Fz  force.  This  is  most  likely 
due  to  the  centrifugal  force  trying  to  pull  the  wing  out  from 
the  sensor.  The  results  from  wing  A  show  a  slight  negative 
Fz  force.  We  conclude  that  the  possible  explaination  for 
the  aerodynamic  force  in  the  negative  Fz  direction  is  wing 
bending  as  shown  in  Fig.  10(A).  The  aerodynamic  force 
acts  normal  to  every  element  of  the  wing  and  it  is  greatest 
close  to  the  wing  tip.  The  resultant  aerodynamic  force 
Ft  bends  the  wing  and  creates  a  negative  Fz  component 
which  seems  to  overcome  the  centrifugal  force  in  this 
particular  case.  This  can  be  further  validated  by  the  fact 
that  the  Fz  component  of  aerodynamic  force  was  found  to 
be  symmetric  between  the  upstroke  and  downstroke. 

Spanwise  twist 

Fig.  9  shows  the  aerodynamic  force  in  the  negative  Fx 
direction.  This  was  found  during  both  the  upstroke  and  the 
downstroke.  The  negative  Fx  force  is  most  likely  due  to  a 
twist  along  the  wing  length  with  the  tip  chord  at  an  angle  of 
roughly  10  to  20  degrees  with  respect  to  the  root  chord.  As 
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Fig.  9.  Aerodynamic  +  inertia  force  of  wing  A,  inertia  force 
for  wing  B  and  filtered  aerodynamic  force  for  wing  A  during 
the  upstroke.  The  unit  of  force  is  Newtons.  A  ten-degree 
polynomial  fit  is  done  on  the  aerodynamic  data. 


shown  in  Fig.  10(B),  the  ith  element  of  the  wing  undergoes 
twisting  by  an  angle  S'fii.  The  aerodynamic  force  dF ^  still 
remains  normal  to  the  ith  element  but  it  contributes  a  force 
in  the  negative  Fx  direction  with  respect  to  the  sensor 
frame. 

A.  Modifications  to  the  aerodynamic  model 

In  order  to  better  match  the  experimental  results,  the 
aerodynamic  model  should  involve  structural  deflections  to 
match  with  the  physical  wing.  As  a  first  approximation,  we 
ignore  the  Fz  component  by  assuming  infinite  rigidity  in 
bending.  We  model  twist  by  assuming  a  linear  variation  of 
twist  from  root  to  the  tip  of  the  wing.  The  twist  5fi>i  at  the 
ith  section  is  given  by 


=  Pji 


0.08 

008’ 


(38) 


where  the  unit  of  r*  and  R  is  meters  and  0.08  is  the 
distance  in  meters  between  the  root  chord  and  the  flapping 
axis  and  /?  is  the  magnitude  of  twist  proportional  to  the 
aerodynamic  moment  about  the  feather  axis.  The  inertial 
moment  is  ignored  due  to  low  flapping  frequency  of  0.5  Hz. 
The  aerodynamic  moment  is  assumed  to  be  proportional  to 
the  normal  velocity  which  is  computed  at  a  point  located 
at  a  distance  of  r  -  along  the  span  from  the  flapping 


Fig.  10.  (A)  Shows  how  the  span  wise  bending  creates  a 

— Fz  component  of  the  aerodynamic  force  which  acts  in  —  Z4 
direction  in  the  sensor  frame.  The  centrifugal  force  acts  in 
the  sensor  +z4  direction.  (B)  shows  how  the  span  wise  twist 
S'lpi  at  the  ith  element  contributes  to  the  —  Fx  component  of 
aerodynamic  force  in  the  sensor  frame. 


axis  and  at  a  distance  of  along  the  chord  from  the 
feather  axis.  Therefore, 

(3  =  ±K0  [un(f2R,t)}2.  (39) 

The  normal  velocity  is  squared  to  give  the  same  effect  as 
aerodynamic  force.  Here,  Kp  is  a  coefficient  to  adjust  the 
magnitude  of  /?  variation  in  a  cycle  due  to  normal  velocity. 
In  the  current  experiment  the  magnitude  of  the  twist  is 
visually  matched  with  the  experiment  by  varying  the  value 
of  Kp.  The  =b  takes  the  sign  of  un(r2R,t)  to  give  the 
correct  direction  of  twist.  The  wing  is  considered  rigid  for 
computation  of  un(r2R,  t). 

Having  an  approximation  of  the  aerodynamic  twist,  the 
total  aerodynamic  force  is  integrated  vectorially  from  the 
root  to  the  tip  of  the  wing.  Considering  the  ith  element  of 
the  wing  located  at  a  distance  ri  from  the  flapping  axis, 
the  total  force  can  be  written  as  follows: 

dFfi  =  Ci^\u(ri,  t)\2 c(ri)dri  + 

C^un{ri,t)c{rifdri.  (40) 

Based  on  the  assumption  that  the  total  force  remains  normal 
to  the  chord  of  any  ith  element,  we  can  transform  this  force 


into  the  sensor  frame  as  follows: 


VI.  Determination  of  Coefficients 


dF*\  (  0  \ 

dFlv  =  /?:  I  dFi  ,  (41) 

dFU  sensor  \  °  J  ^element 

where  R\  is  the  rotation  matrix  between  each  blade  element 
and  the  sensor  frame  given  as 

(C ' S'lpi  0  \ 

Ss^  Csfr  0  .  (42) 

0  0  l) 

Here,  Sipi  is  the  angle  of  twist  of  the  ith  element  with 

respect  to  the  root  chord.  The  moments  in  the  sensor  frame 
can  be  computed  as  follows: 

dMlx  =  -ridFly  ,  dMly=ndFlxy  (43) 


The  sensor  measures  the  total  aerodynamic  force  Ft 
given  by  Eq.  (47).  The  functions  F 2  and  F3  in  Eq. 

(33)  capture  the  physics  of  the  particular  unsteady  effect 
and  the  coefficients  simply  adjust  the  magnitude.  Once  the 
magnitudes  of  C\,  C2  and  C3  are  adjusted  at  any  point 
during  the  flapping  cycle,  these  functions  are  expected  to  be 
robust  enough  to  match  the  experimental  data  throughout 
cycle.  We  determined  the  coefficients  as  follows. 

A.  Determination  of  C\ 

In  order  to  determine  C\,  we  select  the  flapping  kine¬ 
matics  such  that  the  functions  F2  and  F3  become  zero  at 
some  point  in  the  flapping  cycle  but  F\  remains  non-zero. 
This  means  cp  7^  0,  while  ip  =  (p  =  ip  =  0.  The  following 
simple  kinematic  pattern  was  selected  using  Eqs.  (l)-(3). 

(p(t)  =  sin(c of)  ,  (p0  =  (cf>n  =  0,  n  =  1  (48) 


dM\  =  [a  —  d0(ri)\c(ri)dFf ,  (44) 


0(t)=  0  ,  0o  =  0n  =  0  (49) 


where  d0(ri )  is  the  non-dimensional  distance  between  the 
rotation  axis  and  the  wing  leading  edge  as  shown  in  Fig. 
5(A).  The  non-dimenstional  distance  a  =  0.5  gives  the 
mid-chord  location  of  dFf  at  every  section  based  on  our 
assumption  (see  Section  III.  E). 

The  forces  and  moments  can  be  summed  in  the  sensor 
frame  in  order  to  get  the  total  force  and  moment  compo¬ 
nents 


N  N  N 

fx  =  E  dF*’  ^  =  E  dF  ^  =  E dF * :  0  (45) 


N  N  N 

Mx  =  Y^dMl  Mv  =  Y,dMi,  Mz=J2dMi  (46) 

Since  the  Fz  component  of  sensor  force  is  ignored,  the 
resultant  force  Ft  in  the  sensor  frame  can  be  written  as 


ip{t)  =  -4' cos  (ujt)  ,  ^0=  0,  Ci^n  =  T  n  =  1  (5°) 

This  pattern  gives  maximum  values  of  (p  and  ip  at  the 
mid-stroke  position  where  ip  =  0.  The  term  involving  ip 
in  F3  gives  the  rotational  virtual  mass  effect  which  is 
small  in  comparison  to  the  translational  virtual  mass  effect 
and  therefore  we  ignore  it  in  subsequent  analysis.  In  other 
words,  F2  and  F3  both  become  close  to  zero  and  only  the 
translational  force,  Fiev  =  C\F\  is  non-zero  at  the  mid¬ 
stroke.  Therefore,  the  sensor  measures  FT  =  Fiev  and 


Flev 

1/2 p  (p2  r2c(r)dr 


at  (p  —  0, 


(51) 


where  Ft  is  determined  from  force  sensor  components  Fx 
and  Fy.  The  integral  r2c(r)dr  turns  out  to  be  simply 
one  half  of  the  second  moment  of  area  S2  [3].  This  gives 
a  simple  form  for  the  determination  of  C\,  i.e, 


Ft  =  yjFx2  +  Fv2.  (47) 

The  inclusion  of  moderate  spanwise  twist  will  have 
insignificant  effect  on  the  magnitude  of  C\  [5].  We  further 
assume  that  it  has  no  effect  on  coefficients  C2  and  C3. 
Therefore,  based  on  this  assumption,  the  sensor  output,  FT 
could  still  be  used  to  find  the  coefficients  using  Eq.  (33) 
for  rigid  wing  by  ignoring  the  twist  and  Eq.  (40)  and  Eq. 
(45)  can  then  be  used  to  determine  the  force  components 
by  taking  into  account  the  twist.  In  the  computation  of 
kinematics  and  angle  of  attack,  the  twist  Sipi  is  ignored. 


Procedure 


Ci  = 


F*lev 

1/4  p  <P2S2  ’ 


(52) 


For  the  safe  operation  of  robotic  flapper,  the  flapping 
frequency  was  limited  to  0.5  Hz.  This  together  with  the 
wing  geometric  parameters,  gives  the  Reynolds  number  Re 
and  reduced  frequency  K  as  a  function  of  flap  amplitude 


Re  = 


15450tt 


x  <f> 


K  = 


0.2725  x  180  1 


(53) 


180 


7 r 


TABLE  II 
Testing  scheme 


Flap  Amplitude  <E>  (deg). 

Re 

K 

46 

12,404 

0.332 

63 

17,000 

0.2477 

74.5 

20,000 

0.2095 

TABLE  III 

Angle  of  Attack  at  mid-stroke,  (</>  =  0°) 


s.no. 

Rotational  Amplitude  (deg). 

a  (deg) 

* 

90.0 

0.0 

1 

75.0 

15.0 

2 

64.0 

26.0 

3 

52.4 

37.6 

4 

41.2 

48.8 

5 

30.0 

60.0 

6 

19.0 

71.0 

7 

8.0 

82.0 

8 

0.0 

90.0 

Three  different  values  of  flap  amplitude  <f>  were  chosen 
(46,  63,  and  74.5  degrees)  as  shown  in  Table  II.  Each  flap 
amplitude  gives  a  different  value  of  Reynolds  number  Re 
and  reduced  frequency  K.  This  conforms  with  different 
flap  amplitudes  of  flapping  used  by  biological  species. 

For  each  flap  amplitude,  the  rotational  amplitude  4/  and 
eight  correponding  a  from  0  to  90  degrees  were  chosen  at 
the  mid-stroke  position,  as  shown  in  Table  III.  In  fact,  a 
varies  between  -90  and  90  degrees  during  the  cycle,  giving 
two  values  of  force  at  the  same  a ,  differing  only  in  sign. 
We  averaged  the  two  forces  and  determined  an  average 
C i  corresponding  to  a  at  the  mid-stroke  position.  We  did 
not  conduct  the  experiment  at  a  =  0  deg  and  assumed 
Ci(0)  =  0.  This  is  because  normal  force  is  zero  on  a 
symmetric  flat  plate  at  zero  angle  of  attack.  The  chordwise 
friction  force  is  ignored  in  our  model.  In  all,  24  kinematic 
patterns  were  tested. 


Fig.  11.  Coefficient  of  LEV/translational  force  C\  Vs  a 


The  plot  of  Ci  against  a  is  shown  in  Fig.  1 1  for  the  three 
different  stroke  amplitudes.  This  shows  that  the  coefficient 
of  translational  force  varies  linearly  with  a.  Furthermore, 
the  flap  amplitude  and  consequently  Re  has  little  effect  on 
the  slope  (dCi/da).  From  Fig  11,  we  can  approximate  C\ 
as  linearly  varing  with  a  given  by 

Cl  (a)  =  —a,  (54) 

IT 

where  a  is  in  radians. 

B.  Determination  of  C2 

C\  and  C2  occur  implicitly  in  F2.  Therefore,  the  best 
way  to  determine  C2  is  to  adjust  it  untill  the  model  matches 
with  the  experimental  results.  Furthermore,  C2  =  di~d0(r ) 
varies  along  the  span.  dQ(r )  is  known  from  the  wing 
geometry  and  given  in  Table  I.  The  value  of  C2  is  based 
on  di.  We  found  that  di  =  0.75  gives  best  results.  In 
computing  F2,  the  coefficient  C 1  is  considered  a  known 
parameter.  Therefore,  C\  is  determined  before  C2. 

C.  Determination  of  C3 

If  we  modify  the  kinematic  pattern  given  by  Eqs.  (48)- 
(50)  by  taking  4/  =  0°  in  Eq.  (50),  then  F2  becomes  zero 
for  the  entire  cycle  while  F\  becomes  zero  at  the  ends  of 
the  stroke,  i.e,  at  f  =  fmax  and  f  =  0.  However,  F3  is 
maximum  there  since  <fi  is  maximum.  Therefore,  the  total 
force,  Ft,  measured  by  the  sensor  at  the  ends  of  the  stroke 
is  due  to  virtual  mass  effect. 

Ft  F'VirtUal  mass  111  /rr\ 

C3=  -—  =  - - -  ,  at  0  =  (pmax •  (55) 

-G3  -G3 

The  value  of  C3  was  found  to  vary  between  0.5  and  1.0 
for  all  24  kinematic  patterns. 

VII.  Comparison  of  experimental  results 

Most  of  the  published  data  on  flapping  wings  is  for 
flapping  translational  motion  of  the  wing,  i.e;  without 
rotation.  The  only  way  to  compare  it  with  our  experimental 
data  is  to  compare  it  at  the  mid- stroke  of  flapping  cycle 
where  wing  rotation  becomes  zero  instantaneously,  for  all 
the  kinematic  patterns  given  by  Eq.  (48-50).  Secondly,  the 
published  data  is  in  the  form  of  either  coefficients  of  lift 
and  drag  or  coefficients  of  vertical  and  horizontal  force. 
Therefore,  in  order  to  compare  our  experimental  data,  we 
determined  the  coefficients  of  lift  and  drag  at  the  mid- stroke 
position  using  the  lift  and  drag  forces. 

A.  Lift,  drag,  vertical  and  horizontal  force 

The  aerodynamic  force  components  (FXJ  Fyj  Fz)  which 
correspond  to  (£4, 1/4,  £4)  directions  in  the  sensor  frame 
/4,  were  transformed  into  lift,  drag,  vertical  and  horizontal 
components  for  the  entire  cycle  using  Eqs.  (25)  and  (26). 

As  shown  in  Fig.  12,  the  lift  and  vertical  components 
are  positive  for  both  the  upstroke  and  downstroke  while 


Fig.  12.  Lift,  Drag,  horizontal  and  vertical  force  components 
for  <I>  =  46°  and  ^  =  52.4°  in  the  sensor  frame 


the  drag  and  vertical  force  components  cancel  out  during 
the  cycle.  This  is  expected  since  the  kinematic  pattern  given 
by  Eq.  (48-50)  is  symmetric  between  the  upstroke  and 
downstroke.  Therefore,  a  net  vertical  force  is  generated 
corresponding  to  a  hovering  flight.  A  net  horizontal  and 
drag  component  is  produced  if  the  kinematic  pattern  is  not 
symmetric. 

B.  Coefficients  of  lift  and  drag  force 

The  coefficient  of  lift  and  drag  are  computed  from  the 
lift  and  drag  force  plots  at  the  mid-stroke  position  using 
Eq.  (52). 


CL 


Lift 

1/4  p  <j>2S2 


Cd 


Drag 

1/4  p  f2S2 


(56) 


Eq.  (52)  is  applicable  since  Cl  and  Cd  are  the  decom¬ 
position  of  Ci.  The  polar  plot  of  Cl  and  Cd  for  the  three 
flap  amplitudes  (<f>  =  46,  63,  74.5  deg)  is  shown  in  Fig.  13. 
The  polar  plot  shows  high  values  of  Cl  and  Cd  compared 
to  linearly  translating  wings  [5]  and  compares  well  with 
the  published  data  Sane  [11],  Usherwood  and  Ellington  [5], 
[6].  The  only  major  difference  can  be  seen  at  the  maximum 
value  of  Cd  in  Fig.  13  where  the  plot  does  not  go  to  zero. 
This  is  because  we  tested  a  wing  which  could  also  twist 
under  the  aerodynamic  loads.  Therefore,  at  the  maximum 
angle  of  attack,  i.e,  90°  when  the  coefficient  of  drag  is 
maximum  (close  to  3.5),  we  get  added  lift  due  to  span  wise 
twist.  The  polar  plot  also  show  that  the  coefficients  do 
not  vary  with  flapping  amplitudes.  However,  for  a  given 
flapping  frequency,  large  stroke  amplitudes  generate  greater 
lift  due  to  the  fact  that  the  wings  sweep  a  larger  area  with 
higher  translational  velocity.  Therefore,  we  conclude  that 
large  stroke  amplitude  is  vital  for  generating  higher  vertical 
force  during  hovering  flight. 


Fig.  13.  Lift  and  drag  polar  plot  for  three  stroke  amplitudes, 
i.e,  <f>  =  46°, 63°, 74.5°.  These  coefficients  are  computed  at 
the  mid-stroke  position  during  the  flapping  cycle  similar  to 
the  determination  of  C\ 


VIII.  Comparison  of  experimental  results  with 

AERODYNAMIC  MODEL 

The  coefficients  C\,  C2  and  C3  will  now  be  used  in  the 
aerodynamic  model  and  comparison  will  be  made  with  the 
experimental  data.  Note  that  C\  and  C3  were  computed 
from  the  experimental  data  but  C2  was  tuned  to  fit  the 
experimental  results.  In  order  to  compare  the  aerodynamic 
model  with  the  experiment,  we  will  compare  two  kinematic 
patterns.  From  Table  III,  we  select  two  patterns  having 
rotational  amplitude  of  =  75°  and  =  0°  (entries  1  and 
8  in  Table  III).  In  both  cases,  the  flap  amplitude  is  $  =  46°. 
These  represent  the  extreme  cases.  The  experiment  was 
not  conducted  at  ^  =  90°.  If  the  model  compares  well 
with  the  experimental  data  for  these  two  extreme  cases, 
we  will  have  more  confidence  that  the  comparison  will  be 
good  for  the  patterns  in  between.  The  kinematic  patterns 
for  these  two  cases  are  shown  in  Fig.  14,  where  arrows 
indicate  instantaneous  direction  of  motion. 

A.  Comparison  of  kinematics 

In  Fig.  15,  we  compare  the  flap  angle  <f>  and  rotational 
angle  with  the  experimental  data  from  the  encoder. 
The  kinematic  pattern  is  given  by  Eq.  (48)-(50).  The 
complete  cycle  takes  2  seconds  at  flapping  frequency  of 
0.5  cycles/sec.  Note  that  the  maximum  value  of  ^ 
occurs  at  the  mid-stroke  (at  time  =  50  and  150  millisecs) 
with  a  corresponding  minimum  angle  of  attack  during  the 
entire  cycle.  The  comparison  of  <f>  and  is  good  throughout 
the  cycle. 

B.  Comparison  of  twist  5ip 

Fig.  15  also  shows  the  variation  of  twist  5f>  at  the  tip 
section  during  the  cycle  modeled  by  Eqs.  (38)  and  (39). 


Fig.  14.  (A)  This  kinematic  pattern  has  flapping  and  rotational 
amplitudes  of  <1>  =  46°  and  ^  =  74.9°  respectively.  The 
arrows  indicates  instantaneous  direction  of  motion  of  a  wing 
section.  At  mid-stroke,  -0  =  0,  the  section  is  undergoing 
flapping  translation.  (B)  This  kinematic  pattern  has  flapping 
and  rotational  amplitudes  of  4>  =  46°  and  =  0°  respec¬ 
tively.  Section  is  undergoing  translational  motion  only.  Figure 
also  shows  roughly  the  location  of  maximum  twist  Sip  max  as 
observed  during  the  experiment. 


The  model  predicts  maximum  twist  Sip  max  at  the  mid¬ 
stroke  position  (time  =  50,  150  millisecs)  for  pattern  (B). 
This  location  is  also  shown  in  Fig.  14(B).  For  pattern  (A), 
b^max  occurs  when  the  wing  passes  through  the  mid-stroke 
position  and  roughly  close  to  the  end  of  the  stroke  as  shown 
in  Fig  15.  This  location  is  roughly  shown  in  Fig  14(A). 
This  is  because  Sip  is  directly  proportional  to  the  normal 
velocity  according  to  Eq.  (38)  and  (39)  which  turns  out  to 
be  the  greatest  when  the  wing  undergoes  both  rotational 
and  translational  motion  near  the  end  of  the  stroke.  Note 
that  Eq.  (38)  and  (39)  give  an  approximate  representation  of 
the  actual  wing  twist.  The  amplitude  of  twist  was  matched 
visually  with  the  physical  wing  through  trial  and  error  by 
varying  Kp.  During  the  actual  experiment,  we  did  notice 
the  pattern  of  dip  at  the  wing  tip  similar  to  Fig.  15  with 
the  maximum  twist  occuring  near  the  end  of  the  stroke  for 
pattern  (A)  and  at  the  mid-stroke  position  for  pattern  (B). 


Fig.  15.  Comparison  of  kinematic  patterns,  (E)=Experiment, 
(M)=  Model.  Pattern  (A)  <f>  =  46°,  =  74.9°  and  Pattern 

(B)  4>  =  46°,  =  0°  show  a  good  match  with  the  model. 

Plot  also  shows  the  twist  Sip  at  the  tip  section  computed  using 
Eqs.  (38)  and  (39) 


C.  Comparison  of  Fy  component  of  aerodynamic  force 

Fig.  16  shows  the  Fy  component  of  aerodynamic  force  in 
the  sensor  frame  for  the  two  kinematic  patterns.  Plots  (A) 
and  (B)  show  the  force  data  for  the  patterns  having  flapping 
and  rotational  amplitudes  of  <f>  =  46°  and  =  74.9°  while 
(C)  and  (D)  represent  the  pattern  having  amplitudes  of  4>  = 
46°  and  =  0°.  The  Fy  component  was  computed  using 
Eq.  (45).  Figs  16(A)  and  (C)  are  the  plots  of  individual 
aerodynamic  mechanisms  (C\Fi,  F2,  C3F3)  which  are 
components  of  the  total  force  transformed  in  the  sensor  Fy 
direction  for  half  cycle.  Figs  16(B)  and  (D)  shows  how  the 
individual  aerodynamic  mechanisms  contribute  to  the  total 
force  and  match  very  well  with  the  experimental  results. 


Fig.  16.  Comparison  of  Fy  component  of  total  aerodynamic 
force  Ft  —  C\  F\  +  F2(Ci  ,62)  +  C3F3.  (A)  and  (B)  shows 
the  comparison  for  the  pattern  having  amplitudes  <t>  =  46°, 
\p  =  74.9°,  (C)  and  (D)  shows  the  comparison  for  the  pattern 
having  amplitudes  <£>  —  46°,  ^  =  0°. 

For  the  patterns  in  Figs  16(A)  and  (B),  the  rotational 
force  F2  dominates  as  far  as  magnitude  is  concerned.  This 
shows  that  the  quasi- steady  aerodynamic  force  given  by 
C\F\,  is  inadequate  for  modeling  flapping  wing  aerody¬ 
namics  when  the  rotation  is  present.  Furthermore  rotational 
and  virtual  mass  force  F2  and  C3F3  respectively  are  seen 
to  be  zero  at  the  mid-stroke  position  while  the  translational 
force  C\F\  is  zero  at  the  ends.  This  is  because  the  functions 
F2  and  F3  are  zero  while  F\  is  not  zero  at  the  mid- stroke 
(</>  =  0°). 

For  the  pattern  shown  in  (C)  and  (D),  the  rotational 
force  is  zero  throughout  the  cycle  since  ip  =  0  and  the 
wing  undergoes  flapping  translation  only.  Therefore,  we 
see  that  the  trasnlational  force  C\F\  dominates  throughout 
the  cycle  except  at  the  ends  of  the  stroke  where  it  is  zero. 
The  virtual  mass  force  C3F3  is  maximum  at  the  ends  of  the 
stroke  and  adjusts  the  total  force  at  the  ends  as  shown  in 
Fig.  16(D).  This  kinematic  pattern  was  specially  selected 
to  study  the  virtual  mass  effect  since  this  is  the  only  effect 
present  at  the  ends  of  the  stroke. 
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Fig.  17.  Comparison  of  Fx  component  of  total  aerodynamic 
force  Ft .  (A)  shows  the  comparison  for  the  pattern  having 
amplitudes  <E>  =  46°,  =»  74.9°,  (B)  shows  the  comparison 

for  the  pattern  having  amplitudes  4>  =  46°,  ^  =  0° 


D.  Comparison  of  Fx  component  of  aerodynamic  force 

The  Fx  component  of  aerodynamic  force  is  computed 
using  Eq.  (45)  and  compared  with  the  experimental  data  as 
shown  in  Fig.  17  for  the  two  kinematic  patterns.  The  match 
with  the  experimental  data  is  less  accurate  in  the  case  of 
Fx  component.  This  is  probably  due  to  the  approximate 
match  of  twist  5  ft  with  the  actual  wing  in  terms  of  both 
the  spanwise  distribution  as  well  as  variation  during  the 
cycle.  The  model  assumes  linear  twist  from  the  root  to 
the  tip  of  the  wing.  The  actual  twist  distribution  could  be 
different.  It  should  be  noted  that  the  approximation  of  twist 
does  not  affect  the  Fy  component  because  its  magnitude  is 
more  than  three  times  that  of  the  Fx  component  and  the 
error  is  not  visible  in  the  Fy  plot.  We  found  this  typical 
behavior  in  all  the  24  kinematic  patterns. 

E.  Comparison  of  aerodynamic  moments 

The  sensor  measures  the  aerodynamic  moments  as  well 
as  the  aerodynamic  force.  Therefore,  we  can  compare  the 
experimentally  obtained  moments  with  the  moments  com¬ 
puted  from  the  aerodynamic  model.  In  Fig.  18,  the  three 
plots  in  column  (A)  and  (B)  show  the  three  components 
of  moments  ( Mx ,  My,  Mz)  belonging  to  the  two  selected 
kinematic  patterns. 

The  moments  are  computed  using  Eq.  (46).  The  Mx  and 
Mz  components  show  a  good  match  with  the  experimental 
data.  The  accuracy  of  Mx  component  indicates  that  the 
model  accurately  distributes  the  Fy  component  along  the 


Fig.  18.  Comparison  of  moments  ( Mx ,  My,  Mz).  These  are 
computed  based  on  the  assumption  that  the  aerodynamic  force 
at  each  section  dFr  acts  at  the  mid-chord  location.  (A)  shows 
the  results  for  the  kinematic  pattern  having  amplitudes  4>  = 
46°,  =  74.9°,  and  (B)  shows  the  results  for  the  kinematic 

pattern  having  amplitudes  <f>  =  46°,  ^  =  0°. 


span.  The  magnitude  of  Fy  has  already  been  verified  from 
comparison  of  the  force  plots.  Similarly,  the  accuracy 
of  Mz  component  idicates  that  our  assumption  of  mid¬ 
chord  location  of  force  dFT  at  every  section  given  by  the 
parameter  4 a  =  0.5’  in  Eq.  (44)  is  valid.  The  comparison 
of  My  component  is  approximate  due  to  the  fact  that  it 
is  dependent  on  Fx  component  of  force  which  compared 
approximately  with  the  experiment  as  discussed  earlier. 
This  behavior  of  moments  was  found  to  be  consistent  in 
all  24  kinematic  patterns. 

F.  Comparison  of  the  location  of  resultant  aerodynamic 
force 

The  location  of  the  resultant  aerodynamic  force  on  the 
wing  for  one  stroke  is  also  compared  with  the  experiment 
and  shown  in  Fig.  19.  The  distance  a  (chord- wise  direction) 
and  the  distance  b  (spanwise  direction)  give  the  location 
of  the  resultant  aerodynamic  force.  This  is  another  way 
to  test  our  assumption  that  the  total  force  dFT  at  every 
section  acts  at  the  mid-chord  location,  i.e,  a  =  0.5.  a  and 
b  were  calculated  using  the  values  of  force  and  moment  at 
the  location  of  force/torque  sensor  for  both  the  model  and 
experiment  using  the  following  equations. 


(57) 


Fig.  19  shows  good  comparison  of  aerodynamic  and  exper¬ 
imental  values  of  force  location  (a  and  b)  for  one  kinematic 
pattern.  Fig.  19  also  indicates  that  the  location  of  force  does 
not  vary  during  the  stroke  except  near  the  very  ends. 


Fig.  19.  Location  of  resultant  aerodynamic  force  on  the  wing 
given  by  the  distance  a  and  b  for  the  kinematic  pattern  (<f>  = 
46°  and  A/  =  64°).  Note  ‘Exp’  stands  for  experimental  and 
‘Mdl’  stands  for  aerodynamic  model.) 


IX.  Conclusion 

This  paper  presents  a  methodology  for  the  experimental 
determination  of  unsteady  aerodynamic  force  coefficients 
based  on  the  principle  of  dynamic  similarity.  These  coef¬ 
ficients  are  used  in  the  quasi-unsteady  aerodynamic  model 
which  additionally  takes  into  account  the  wing  twist  due  to 
aerodynamic  loads.  The  key  conclusions  drawn  from  this 
work  are  given  below 

1)  We  found  a  good  match  of  Fy,  Fz,  Mx  and  Mz 
components  and  an  approximate  match  of  Fx  and 
My  components  of  resultant  force  and  moment  at  the 
base  of  the  wing  for  all  the  24  kinematic  patterns. 

2)  The  coefficient  of  translational  force  C\  is  deter¬ 
mined  as  a  function  of  angle  of  attack  at  one  point 
during  the  flap  cycle.  Similarly,  the  coefficient  C2 
and  Cs  are  determined  as  constants  for  the  entire 
flapping  cycle.  However,  the  model  matches  the 
experimental  data  for  the  entire  cycle  and  for  all 
kinematic  patterns.  This  idicates  that  the  functions 
Fi,  F2  and  F3  are  robust  and  capture  the  trend  or 
physics  of  the  unsteady  effects. 

3)  The  variation  of  stroke  amplitude  <f>  at  constant 
flapping  frequency  does  not  alter  the  translational 
force  coefficient  C\.  However,  an  increase  in  stroke 
amplitude  results  in  higher  lift  and  drag  forces  during 
the  cycle.  This  means  that  high  stroke  amplitude  is 
essential  for  hovering  flight. 

4)  The  model  assumes  that  for  high  angles  of  attack, 
typically  encountered  in  hovering  flapping  flight,  the 
total  aerodynamic  force  remains  normal  to  the  chord 
and  the  chordwise  component  is  zero  for  the  entire 
cycle.  The  good  match  of  model  and  experimental 
results  indicates  that  the  total  force  remains  roughly 
normal  to  the  chord  for  the  entire  cycle  as  assumed. 
In  making  this  conclusion,  we  have  taken  into  ac¬ 
count  the  approximate  span  wise  twist  in  the  model. 


5)  We  also  found  that  the  effects  of  virtual  mass  force 
are  small  compared  to  translational  and  rotational 
effects  at  the  Re  range  of  12,000  to  20,000.  In  our 
model,  we  have  assumed  that  C3F3  is  virtual  mass 
effect.  However,  it  might  be  wake  capture  effect 
because  both  act  in  a  similar  manner.  This  needs  to 
be  verified  using  PIV  measurements. 

6)  Moderate  spanwise  twist  has  very  little  effect  on  the 
magnitude  of  C\  and  total  force  Ft ,  however,  it 
changes  the  force  components  in  the  sensor  frame 
significantly.  This  implies  significant  change  in  the 
direction  of  resultant  force  vector.  For  wings  with 
large  spanwise  twist,  the  procedure  should  be  mod¬ 
ified  by  decomposing  C\  into  horizontal  Cp  and 
Vertical  Cv  coefficient  of  force.  However,  we  have 
not  pursued  this  method. 
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